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ABSTRACT 

It has been known since the 1950's that the observed gas content of Galactic globular clusters 
(GCs) is 2 — 3 orders of magnitude less than the mass lost by stars between Galactic disk 
crossings. In this work we address the question: What happens to this stellar gas? 

Using an Eulerian nested grid code, we present 3D simulations to determine how stellar 
wind material evolves within the GC environment. We expand upon work done in the 70's and 
move a single-mass King-model GC through the Galactic halo medium, stripping a 10 5 M Q 
GC of its intra-cluster medium but predicting a detectable medium for a 10 6 M cluster. We 
find from new multi-mass King model simulations, the first to incorporate empirical mass- 
loss formulae, that the single-mass King model underestimates the retention of intra-cluster 
gas in the cluster. Lastly, we present a simple discretised multi-mass GC model, which yields 
lower levels of intra-cluster medium compared to the continuous single- and multi-mass King 
models. 

Our results show that there is still an issue with the predicted intra-cluster gas content of 
massive GCs. We conclude that by modelling GC systems more accurately, in particular the 
stellar structure and description of mass loss, we will be able to work towards resolving this 
issue and begin to fill in some of the gaps in our understanding of the evolution of globular 
clusters. 

Key words: Galaxy: globular clusters: general - stars: mass-loss - stars: Population II - ISM: 
evolution - hydrodynamics - methods: numerical. 



1 INTRODUCTION 

This paper addresses the long-standing problem of why we observe 
no, or very little, intra-cluster medium (ICM) gas in Galactic glob- 
ular clusters, when it is understood that giant branch stars are con- 
tinuously supplying enough material to form an ICM. 

Galactic globular clusters (GCs) are massive star clusters 
(10 4 — 1O 6 M0) of approximately 13Gyr in age and are dis- 
tributed roughly spherically about the Milky Way centre with dis- 
tances ranging from several to a few tens of kpc and some as far 
as 100 kpc. They orbit the galaxy on timescales of a few 10 s yr 
(Odenkirchen et al. 1997) and will cross the Galactic disk twice 
per orbit, removing ICM material from the cluster potential each 
time. The more distant clusters have orbits of a few Gyr with 
some perhaps never crossing the DislQ There is a wealth of obser- 
vational evidence for post-main sequence mass loss in GCs from 
dust tracers (|Origlia et al.|1997||R amdani & Jorissen 2001 ; Origlia 
et al.|2002| McDonal d et al.|2009} and gas tracers (|Cohen|1976| 
Dupree|1986||Bates et al.|1993||McDonald & van Loon|2007V The 

* e-mail: wrp@astro.livjm.ac.uk 

1 The Disk and Halo, capitalised, refer to components of the Milky Way. 



winds from these evolved giant branch stars will supply material to 
the GC ICM on timescales much less than the GC orbital period. 
Over the whole red giant branch (RGB) phase, stars typically lose 
0.2 — 0.3 Mq of gas in the form of cool, slow winds, with veloc- 
ities around 10 — 20 kms -1 , which is below the escape velocities 
for GCs (25 — 75kms~ 1 ). Asymptotic giant branch (AGB) stars 
lose yet another ~ 0.1 Mq with velocities of the same order as the 
RGB winds. Depending on the mass of the GC, it can be expected 
that 10 - 100 M are available to the GC ICM between Galactic 
disk crossing events (Tayler & Woodj, 1975| >. 

1.1 Discrepancy between naive expectations and 
observations 

Since the late 1950's, there have been ongoing searches for an ICM 
that resides within GCs; they have mostly focussed on detecting 
hydrogen (molecular, atomic, ionised) or dust (see Roberts 1988; 
|van Loon et al.|2006| for a summary). These observational studies 
show two orders of magnitude less gas than expected in any de- 
tectable form. Most attempts only place upper limits on the mass 
of ICM gas with one or two yielding tentative detections. Here, we 
provide a brief review of the approaches taken towards detecting an 
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ICM. 

Studies targeting molecular hydrogen (via CO, OH and H2O) 
yield upper limits of the order of 0.1 Mq or no detection fKnapp 
& Kerr|[l973l |Kerr et al.|[T976l |Cohen & Malkanl[T979l |Frail &| 



Beasle yfT994l|Smith et 5JT9 95 ; Leon & Combes|1996||van Loon| 
et al.|2006) . Efforts focussed on atomic hydrogen using the 21 cm 



line give upper limits from a few to a few tenths of a solar mass 
l|Heiles & Heirry|[T966l |Kerr & Knapppln^ [Knapp et al.||1973| 
Conklin & Kimble|1974||Kerr et al.l 19761 ISniith et al.|1990l|van 



Loon et al.|2006). The most sensitive observation to date has been 



van Loon et al. ( 200 9}, in which no gas is detected in the 



that of 

four clusters they studied. |van Loon et al. 12006) give a tentative 
detection of M H ~ 0.3 M in M15. 

Ionised hydrogen searches (via Ha and continuum free-free 
emission) claim upper limits in the range ~ 0.1 — 1 Mq jSmith| 
|et al.|1976||Faulkner & Freeman|1977||Knapp et al.|1996} . A com- 
plementary studybylFreire et al. (2001} infers a population of free 
electrons in the cluster 47 Tucanae from pulsar timing observations, 
indicating the existence of 0.1 Mq of ionised gas within that clus- 
ter. 

In addition to investigations concentrating on hydrogen, there 
have also been sub-millimetre and infra-red observations focussing 
on dust ( Knapp e t al.| 1 995 1 |Hopwood et al. 1 1 998]| 1 999| [MaT sunaga 
|et al.|2008||Boyer et al.|2008||Barmby et al.[2009} . These give up- 
per limits on the dust mass of M duat ~ 10 - 10~ 5 M Q . M15 is 
the only cluster with a detection of dust, where Evans et aTT] ( 2003 I 
determine M duBt ~ 5 x 10~ 4 M and |Boyer et al. 1(2006} find 
M dust ~9 x KT 4 M . 

Clearly, some mechanism is at work for removing ICM gas 
from the GC environment. 



1.2 Two theoretical avenues for resolution 

Theoretical investigations into removing gas from GCs can be 
roughly split into two approaches: mechanisms intrinsic to the GC 
and those arising from the GC environment. The following is an 
overview of the theoretical picture to date. 

Initial investigations looking at intrinsic mechanisms solved 
steady-state flow equations in one dimension (Burke 1968; [Scott] 
& Rose 1975; Faulkner & Freeman 1977). Faulkner & Freeman 



1977) advanced the work of [Burke (1968} and |Scott & Rose 
1975 ) and included energy injection by stellar winds and ra 



diative cooling in the ICM gas. The solution of the steady-state 
equations require larger RGB wind velocities than are still cur- 
rently observed. The conclusions warrant further observational ef- 
forts and provided motivation for time-dependent simulations of 
GC ICM gas. The first hydrodynamical simulations ( Van denberg| 
& Faulkner|l977t traced the evolution of ICM gas in a 10 b M 



GC using RGB wind velocities not permitted by steady-state mod- 
els (wRGBwind < 118 km s -1 ). These lower RGB wind simula- 
tions predict that ICM gas will sink to the centre of the cluster 
and cool. Vandenberg ( 1978 ) included heating of the ICM due to 
a UV radiation field from hot horizontal branch stars. This lowered 
the steady-state wind model threshold RGB wind velocity from 
WRGBwind ^ 118 kms -1 to tiRGBwind ^ QS kms" 1 , aga i n pro- 
viding support for further observational efforts. Vandenberg ( 1978 ) 
notes that not all GCs possess hot HB stars. Coleman & Worden 
{1977} included the flaring activity of M-dwarf stars to produce an 
outflow of ICM gas. The caveat however, is that the flaring proper- 
ties, numbers and distribution of M-dwarfs within GCs are highly 
uncertain. Furthermore, modelling the M-dwarf flaring interaction 
with ICM gas requires additional assumptions. Scot t & Durisen| 



(1978) produced steady-state models of ICM outflows driven by 
novae, assuming radiative cooling to be unimportant. However, 
there are several uncertain factors and novae in GCs are much 
less common than previ ously expected (Dobrotka et al. 2006| |Ma7>| 
|carone & Knigge|2007| with an estimated rate of 4 x 10~' i yr~ 
(Bode & Evans 2008). Furthermore, novae have axisymmetric out- 
flows (O'Brien et al. 2006) and will be less efficient at removing 
gas from the cluster potential. 

Extraneous to simulations, there are several studies that have 
investigated potential gas removal mechanisms that may be active 
within GCs. Spergel ( 1991 1 looked towards energy liberated from 
pulsar winds whilst Yokoo & Fukue (1992 ) did the same for X-ray 
bursters. These letters are brief and only assess the energy require- 
ments for lifting gas from a GC potential well. They contain no 
detailed investigation as to how energy from these objects is trans- 
ferred to the ICM gas. |Umbreit et al^] ( |2008[ l propose stellar colli- 
sions as a method for removing gas in M15 and |Dupreeet~aT] (2009) 
suggest that low-metallicity RGB winds may be intrinsically fast 
enough to leave the cluster potential. 

Extrinsic mechanisms influencing ICM evolution generally 
concern the effect of the GC moving through the Galactic halo 
medium ( prank & Gisler 1976). Such interactions cannot be in- 
cluded in ID simulations or steady-state models. These authors 
used an analytical equation to assess the Galactic halo properties 
required in order to strip a massive GC of ICM gas. They calculated 
a density roughly an order of magnitude higher than that observed 
in the Halo and the halos of other galaxies (Spitzer 1956, Silk 1974; 
Fukugita & Peebles 2006 1. To date, no numerical simulations into 
the Galactic halo-GC ICM interaction have been performed. 

The questions remain: what happens to gas from giant branch 
stars and what causes the low ICM content in GCs despite this mass 
loss? This paper presents the first 3D hydrodynamical simulations 
to trace the evolution of gas within the globular cluster environ- 
ment. We model the GCs motion through the Galactic halo and 
study the effect it has on the ICM gas evolution (which is not pos- 
sible in spherical symmetry). We provide a more realistic descrip- 
tion of the GC stellar population via a multi-mass King model and 
apply empirical mass-loss laws to describe mass loss within the 
cluster (from RGB and AGB stars). Using the multi-mass and King 
models, we also compare the consequence of altering the radial dis- 
tribution of the stellar mass loss within the GC. We further refine 
the multi-mass GC into a discrete population of point masses orbit- 
ing within the GC potential. 

This paper is organised as follows. We justify the methods 
used together with a description of our GC model and present 
the validation of our hydrodynamics code against Vanden berg &| 
|Faulkner|jl977F in Section [2] In Section [5] we discuss the estab- 
lishment of the GCs motion through the Galactic halo medium. 
We justify and describe the use of a multi-mass King model with 
empirical mass-loss laws in Section|4]and continue to examine the 
effect of discretising the multi-mass model in Section [5] We con- 
clude the paper in Section [6] with a brief summary and discussion 
of further work that can be done with this kind of investigation. 



2 NUMERICAL METHODS AND VALIDATION 

2.1 Hydrodynamics, initial conditions and free parameters 

GCs are essentially spherically symmetric systems where the stel- 
lar structure is well represented by ID models. However, at least 
a 2D simulation is required in order to simulate the interaction of 
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the ICM gas with the Galactic halo. For the case of a discrete stel- 
lar population, a fully 3D simulation is essential. We take a 3D 
Eulerian hydrodynamics code (derived from PROMETHEUS) that 



|Vandenberg & Fa ulkner ( 1977 1 where p, v and e are the ICM gas 
density, velocity and specific energy, respectively. U a is the total in- 
jected specific energy and contains only the specific kinetic energy 



ward 



1992 



employs the piecewise parabolic method (PPM) ^Colella & Wood-] injected from the expanding stellar wind, /3 (= §Uwind)> an d the 
1984) and utilises a fixed nested grid architecture (see[Ruffert change in bulk specific kinetic energy as the injected wind comes 

to rest relative to the flow of the ambient ICM (i.e. after momentum 
conservation). VF77 ignore the contribution to U a from the thermal 
energy of the stellar winds and so it is not included in our tests. U a 
is given by 



for details). The density scale of the stellar population in GCs 
spans many orders of magnitude, from the core to the tidal radius. 
Therefore, nested grids allow us to resolve the dense central regions 
and cover the whole cluster with minimal computational cost. 

We follow the same procedure as previous studies in modify- 
ing this code (Faulkner & Freeman 1977, Vandenberg & Faulkner 
1 1977 1. We employ the isotropic (i.e. non-rotating) King model 
(King 19661 in order to represent the stellar structure of a globu- 
lar cluster in our simulation volume. The King model assumes a 
single homogeneous stellar population (i.e. equal-mass stars), de- 
scribed by a finite isothermal sphere. The distribution function (DF) 
is a Maxwellian one minus a constant (Equ. {I); see King 1966, for 
details) 



/I 



(V-Vo) 



(1) 



where k is a dimensional constant, V is the local gravitational po- 
tential (where Vo denotes the central value), a is the stellar velocity 
dispersion and v* is the local stellar velocity (where w» e denotes 
the local escape velocity). The integration of this DF via the Pois- 
son equation produces a stellar density profile that successfully re- 
produces the surface brightness profiles of many Galactic GCs (i.e. 
a flattened core and a cut-off tidal radius where the density falls 
to zero). In addition to the stellar density profile, the integration of 
the DF also yields the radial profiles of the gravitational potential 
and mean squared stellar velocity. The purpose of this stellar dis- 
tribution is to inject gas into the inter-cluster medium. Stellar mass 
loss is typically prescribed via the specific mass-loss rate, a, that is 
treated as a free parameter. The product of a and the local stellar 
density, p„, gives the local stellar mass-loss rate per unit volume, 



at p* 



(2) 



The GC model is placed centrally in the simulation volume and 
implemented by prescribing each cell with a local stellar density, 
gravitational potential and mean squared stellar velocity, according 
to its radial position. The hydrodynamic models are advanced in 
time explicitly and observe the following sequence of steps: the 
stellar wind material is first injected onto the grid, conserving local 
mass, 



momentum, 



Pn+l = Pn + Pinj , 

Pn + l 1>n+l = Pinj (l>* + «wind) + pn V„ 
= pnV n , 



and energy, 



Pn+l Sft 



(3) 



(4) 



(5) 



Pn+l e n +l = e n p n + UaPinj 

accordingly. The momentum term in Equation ^ for the injected 
gas vanishes since we assume a spherically symmetric stellar wind 
and, for an isotropic stellar velocity field, v* = when integrated 
over all velocity vectors for the 'stars' at a given point. After con- 
serving mass, momentum and energy, radiative cooling (5ft) is ap- 
plied before advancing the grid hydrodynamically. The subscripts 
n and n + 1 refer to the current model and the one being advanced 
in time by At, respectively, whilst inj and * refer to injected gas 
and stellar properties, respectively. We adopt the same notation as 



U a 



Vn+1 



(6) 



where the cross term, v» ,v n +i from expanding the quadratic disap- 
pears due to the isotropic stellar velocity field, (vf) is the local (i.e. 
dependent on radial position) mean squared stellar velocity and is 
derived from equation (31) in King ( 19661. 



(vi) 



4 2 
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rW 

Jo e " 



3 , 



(7) 



Equation |7J appears, though with a typographical error, as equa- 
tion (5) in Faulkner & Freeman {1977). W is the dimensionless 
gravitational potential and relates to the gravitational potential via 



W = 



V 



(8) 



and r\ is the dimensionless stellar specific kinetic energy, denoted 
by 



n = 



2<r 2 



(9) 



The denominator in Equation |7| can be integrated by parts to pro- 
duce a power law series in W, 



E 



n 

m— 3 



2m - 1 



(10) 



2 s 22 i 222 a 



We employ t he same King model parameters as Van denberg] 



& Faulkner 



\\91l\ (hereafter VF77) to produce a If? M GC. 
The parameters adopted are: core radius r c = 0.5 pc, tidal radius 
r\ — 23.9pc, central stellar density p = 1.738 x 10~ 17 gcm -3 
and central dimensionless gravitational potential Wo = 7.5. We 
adopt the same metallicity (X = 0.7, Y = 0.298, Z = 0.002) and 
radiative cooling of the ICM as well. The treatment of the radia- 
tive cooling is described in Faulkner & Freeman ( 1977 1 and is an 
analytical equation (their equations (11), (12), (13) and (16)) that 
approximates a cooling function that models free-free emission, 
the collisional excitation of permitted lines as well as forbidden 
and semi- forbidden lines ( |Cox & TuckeT|[T969l |Cox & Daltabuitl 
|1971[ >. The specific mass-loss rate (for RGB winds) takes the value 
a = 4 x 10~ 19 s _1 . VF77 perform three simulations to probe the 
effect of the stellar wind velocity on the ICM evolution in GCs. 
The adopted velocities are 50, 100 and 150 km s" 1 ; faster than 
that expected for RGB winds. Stellar wind velocities are a source 
of energy, injected into the ICM via the specific wind kinetic en- 
ergy, p. Our validation simulations are performed on seven nested 
grids, each a 64 3 Cartesian mesh. The main (i.e. largest) simulation 
volume has a side length that is four times the tidal radius. Seven 
nested grids enables data to be produced at the same length scales 
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Table 1. Summary of results for the 50kms 1 wind simulation at 
0.2 Myr. 



ICM properties at 0.2 Myr 


this work 


VF77 


log (po/ gem -3 ) 


-21.26 


-20.5 


log (To/ K) 


3.99 


3.93 


^radial, max/ km s 


11 


5 


^radial, mill/ km s * 


-22 


-30 


log (f reservoir/ Po) 


-0.175 


-0.51 



in the core as VF77. The initial conditions are a stationary Galac- 
tic halo medium with a uniform density of pn = 10~ 27 gem -3 
at a temperature of Th = 10 5 ' 5 K. We do not attempt to confirm 
VF77's low pressure Halo simulations where Th = 10 2 K. 



2.2 Comparison of our results with those of Vand enberg & 
|Faulkiier|p77) 

We validate our hydrodynamics code with the reproduction of the 
ID simulations of VF77 in 3D. We perform the same three sim- 
ulations for a 10 Mq GC, as presented in VF77, on our simu- 
lation volume. In this section, we split the analysis of the results 
into two discussions: the unsteady cluster wind simulations (i.e. ft 
corresponding to 50kms _1 and 100 km s -1 ) and the steady-state 
cluster wind simulation (J3 analogous to 150 km s -1 ). 

In order to save computing time, VF77 begin their models af- 
ter some time At, performing analytic calculations to determine the 
initial density, temperature and velocity profile for the gas on their 
ID grid. This neglects the effects of cooling in the initial model. 
Since VF77 do not provide the value of At, we start our simula- 
tions from t = 0. 



2.2. 1 Unsteady cluster wind simulations 

VF77's results show that the 50 km s" 1 and 100 km s' 1 RGB wind 
velocities result in gas failing to escape the cluster potential, lead- 
ing to the growth of the ICM within the GC core. They present 
several models from different times in their simulations (see fig. 5 
and fig. 3, respectively, of VF77). A selection of the correspond- 
ing models from our simulations are presented in Fig.[TJand Fig. [2] 
respectively. Due to limited computer resources at the time of vali- 
dation, we only simulate the minimum number of models required 
to perform a detailed comparison with VF77. As can be seen in Fig. 
[T]and Fig. [2] the main morphological features of the VF77 models 
are reproduced by our code. The gas density in the central regions 
increases, forming a bound sphere which cools over time. Gas ex- 
ternal to this central reservoir falls onto the surface of the spherical 
region except at larger radii where gas is seen to move outwards 
from the cluster. Unsurprisingly, a more detailed comparison with 
VF77 yields some differences in the gas properties between equiv- 
alent models. 

For the 50kms _1 simulation, we compare the gas properties 
of some notable features from our model e and VF77's model E 
(their fig. 5) in Table [T] We find our central density to be almost 
a sixth of VF77's and the extent of the central reservoir of gas is 
twice the size. The central temperature is 15% higher, likely the re- 
sult of a lower cooling rate in the lower-density reservoir. The peak 
inflow velocity from VF77, shows a higher inflow velocity due to 
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Figure 1. A selection of unsteady wind models from our simulations of a 
10 6 Mq GC with an assumed f w i n d = 50kms _1 . Model e— >f shows 
a model between those of E and F in VF77's fig. 5. The inner and outer 
dashed vertical lines show the position of the core radius and tidal radius, 
respectively. 



gas falling onto a more compact core. The outflow velocity feature 
coincides at the same radii, yet our results have a velocity that is 
just over twice that of VF77. 

For the 100 km s -1 wind simulations, we compare model e 
with model E of VF77 (their fig. 3). In Table|2] the central density 
is nearly a sixth the value in VF77 and the central temperature is 
about 19% higher. The radial extent of our core reservoir of gas is 
about 33% larger. The peak inflow velocity is half that of VF77, 
and occurs at a slightly larger radius. The outflow peak velocity is 
just 6% less than VF77's. 

VF77 do not publish their ICM gas energy data, therefore, the 
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Table 2. Summary of results for the 100 kms 1 wind simulation at 
0.1 Myr. 
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Figure 2. A selection of unsteady wind models from our simulations of a 
10 6 Mq GC with an assumed u W j n< j = 100 kms -1 . Model e— >f shows 
a model between those of E and F in VF77's fig. 3. The inner and outer 
dashed vertical lines show the position of the core radius and tidal radius, 
respectively. 



temperature data presented in Tables[T]and[2]are for reference pur- 
poses only. We warn the reader that it is inadvisable to directly com- 
pare the gas temperatures because, unlike energy, it is not solved for 
in the conservative form of the hydrodynamics equations. 

Aside from the initial models used in the corresponding sim- 
ulations, there are other systematic differences leading to a diver- 
gence in the solutions. The implicit method employed by VF77 and 
our explicit scheme are fundamentally different and, naturally, will 
not lead to exactly the same numerical solution. Furthermore, a ID 
simulation requires a central boundary condition. VF77 apply the 



ICM properties at 0.1 Myr 


this work 


VF77 


log (Po/ gem -3 ) 


-22.5 


-21.75 


log (To/ K) 


4.1 


4.025 


^radial, max/ ^ m s 


44.75 


47.5 


^radial, min/ k m s 


-8.75 


-17.5 


log ('^reservoir / P^) 


-0.575 


-0.7 



Table 3. Summary of features in the steady-state model at 0.5 Myr for the 
150 kms" 1 wind simulation. 



ICM properties at 0.5 Myr this work VF77 

log (po/gem" 3 ) -23.95 -24.01 
log (To / K) 5.37 5.48 
fradial (rtj/kms" 1 109 117 
Micm/Mq 2.83 2.59 



commonly adopted condition that v — at r — 0, therefore, gas 
approaching the centre must become stationary and cannot flow 
freely around this point. An excellent example of this restriction 
can be seen in model d of Fig. [5] where the gas has a velocity 
of — 2.5 km s -1 at very small radii. In contrast, curve D of fig. 3 
(VF77) shows this negative velocity approaching zero at smaller 
radii. We argue that this boundary condition artificially traps gas 
at smaller radii leading to a more condensed and cooler reservoir 
of ICM gas. With these in mind, and considering that the initial 
density is five to six orders of magnitude lower, the ~ 0.75 dex 
difference in the central density (Tables [T] and |2j is not as alarm- 
ing as it first seems. Considering these systematic differences, it is 
encouraging that the morphology is faithfully reproduced by our 
code. 



2.2.2 Steady-state cluster wind simulation 

The 150 kms -1 RGB wind simulation leads to a steady-state out- 
flow of ICM gas from the globular cluster (see VF77, fig. 1). VF77 
compare their GC wind model (curve E) with the steady-state so- 
lution of Faulkner & Freeman ( 1977 ) (the dotted line in their fig. 
1) and find them to be in excellent agreement. Since the evolution 
of the ICM leads to a steady-state solution, the central boundary 
condition exerts little influence on the gas dynamics. Furthermore, 
cooling does not dominate the flow for this RGB wind velocity, re- 
ducing the number of factors to consider compared to the unsteady 
GC wind case. Consequently, we therefore expect an overall closer 
agreement with VF77. 

Fig. [3] shows models from our simulation and is the equiva- 
lent of fig 1. in VF77. Since a GC wind is established within rela- 
tively short timescales, we were able to reproduce the same models 
presented in VF77. Encouragingly, the two sets of models are re- 
markably similar. However, the temperature profile of models a, 
b and C, show a striking disparity, which is less so in later mod- 
els. Our models show a lower temperature, at low and interme- 
diate radii, compared to VF77. This difference, is caused by the 
absence of cooling in VF77's method for determining their initial 



6 W.R. Priestley, M. Ruffert, M. Salaris 




80 - 



(3^)!/= - 150 kms"' 
Model Time (yr) 
a 5.00 x 10' 
5.22 x 10 3 
2.08 x 10 4 
1.00 x 10 5 
5.00 x 10 6 







hlM 




MM 


M M M M 1 




M M M 




-1 1 

log(r/pc) 



Figure 3. A selection of the steady-state wind models from our simulations 
of a 10 6 Mq GC with an assumed f w i n( j = 150 kms" 1 . The inner and 
outer dashed vertical lines show the position of the core radius and tidal 
radius, respectively. 



model. To illustrate this, we put AT = 50 yr (i.e. the time to reach 
model a) into equations (14) and (15) of VF77. The corresponding 
central density and temperature are log (p / gem" 3 ) = —25.92 
and log (To/ K) = 5.73. The estimated density agrees, to within 
5%, with our model a, however, the analytical temperature is 58% 
higher than ours. Compared to VF77's central gas temperature of 
log (To / K) = 5.69, our model a has a temperature that is 31% 
lower at log(T /K) = 5.53. Table [| gives a quantitative com- 
parison of some notable features in model e, the model showing 
a steady-state outflow of gas from the GC. The central density we 
obtain is about 15% higher than VF77, and the central temperature 



22% lower. At the tidal radius, our radial gas velocity is about 6% 
lower than VF77. Due to the denser, cooler central gas parameters, 
our ICM is 9% more massive. Despite some differences in the early 
models, the overall morphology at steady state is almost identical, 
with a quantitative agreement to within 10 — 20%. 

2.2.3 Summary of comparison 

In summary, we run simulations employing the same GC models, 
gas physics, Halo gas properties and RGB stellar wind velocities, 
as presented in Vandenberg & Faulkner ( 1977 ). We compare the re- 
sults in order to validate our hydrodynamics code. In all cases, the 
ICM evolution is qualitatively reproduced, with equivalent models 
showing few morphological differences. We identify several sys- 
tematic quantitative differences in the results, which we attribute to 
several factors: the fundamental difference between an implicit and 
an explicit code in reaching a solution; a ID mesh requires a central 
boundary condition of v = at r = that will exert a "braking" 
influence on the gas dynamics in the central regions; VF77's ini- 
tial models start after some undefined time period At and ours at 
t = 0. In spite of these points, our code reproduces the ICM evolu- 
tion presented in VF77 remarkably well. 

To a crude approximation, simulating from t — can be 
taken to represent the environment shortly after a catastrophic ICM- 
removal event. 



3 MOVING THROUGH THE HALO ENVIRONMENT 

The GC model in section [2T2] was at rest with respect to the ambi- 
ent medium. Obviously GCs are not stationary within the Galactic 
halo but in continuous motion throug h this hot, ten uous ambient 
gas (T ~ 10 6 K, p ~ 10~ 27 gcm~ 3 ; |Spitzer|l956|> at supersonic 
speeds between 100 kms" 1 and 300 kms - lOdenkirchen et al. 



|1997| >. |Frank & Gisler| ( |1976) writ e that their analytical formula is 
confirmed by the results of |Gisler| ( fl976| >, who undertakes 2D sim- 
ulations of elliptical galaxies in clusters. There have been no direct 
simulations investigating a GCs motion within the Galactic halo in 
the literature to date. In order to determine if such simulations are 
justified, we use equation (2) from Frank & Gisler ( 19761 to work 
out the Halo density required to continuously strip our model GC 
of its ICM gas. We roughly estimate the surface density of stars at 
the centre of the cluster, er , to be p r c = 134kgm" 2 . We use 
vqc = 200 kms - as a typical GC velocity, which is the me- 
dianof the distribution of values presented in Odenkirchen et al7] 
1 1997 i. The escape velocity for our clust er is v esc — 76kms - 
and q = 4 x 10 -19 s -1 , as in Section 2.1 The corresponding Halo 
density, must be around 10 -25 gem" 3 in order to remove all gas 
from the GC; two orders of magnitude higher than the upper Halo 
value given by Spitzer ( 1956). Therefore, away from the Disk, we 
should still expect GCs to build-up some form of ICM. 



3.1 What environmental factors are important? 

We perform simulations of the GC moving through the ambient 
Halo medium. This is achieved by following the GCs motion in 
the GCs frame of reference and choosing a set of boundary con- 
ditions that produces a continuous flow of ambient Halo gas past 
the GC. The GC remains fixed at the centre of the simulation vol- 
ume. In all simulations, gas flows in the positive x direction of the 
grid. We run several low-resolution simulations, using three nested 
grids at a resolution of 64 3 , in order to test the influence of several 
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Table 4. Initial conditions for the low-resolution simulations of a GC mov- 
ing through the Galactic halo medium. 



Simulation 


M G c 


PH 


fGC 


a 




(Ms) 


(gem" 3 ) 


(kms" 1 ) 


(s-i) 


A 


10 s 


10 -26 


200 


4 x IO" 19 


B 


10 6 


10 -26 


200 


4 x 10~ 20 


C 


10 B 


10 -26 


200 


4 x IO" 19 


D 


10 6 


io- 27 


200 


4 x 10~ 19 


E 


10 6 


io~ 27 


200 


4 x 10" 20 


F 


10 s 


ur 27 


200 


4 x 10~ 19 


G 


10 6 


IO" 27 


300 


4 x 10~ 19 



parameters on ICM evolution. These parameters are, the GC veloc- 
ity, the specific mass-loss rate (a) and the density of the Galactic 
halo medium. A convergence test showed, that using three nested 
grids rather than seven still gives the same overall results regarding 
ICM gas physics and morphology, but slightly overestimates the 
ICM gas content. However, an understanding of the influence of 
different parameters over longer evolutionary timescales is worth 
the trade-off in accuracy of the predicted ICM mass. The GC ve- 
locity is set to either a median 200 km s -1 or an upper value of 
300 kms -1 . The specific mass-loss rate, a, takes eith er one of the 
two extreme values 4 x 1 0~ 19 s~ 1 (determined in Knapp et al. 

s" 1 jTayler & Wood|l975| >. Finally, in ad- 



1973| and 4 x IO -20 

dition to the assumed 10 _/ ' gem -3 , we also give the Halo den 
sity an extreme value of 10 -26 gcm -3 in some simulations. We 
provide a list of the initial conditions used for each low-resolution 
simulation in Table [4] The stellar wind velocity is assumed to be 
20 kms -1 for all simulations. 



3.2 ICM-Halo gas interaction 

We ran simulations for at least lOMyr (i.e. 10% of the time be- 
tween Disk crossings). We also compare models from different 
simulations at this time. We first compare the ICM gas mass evolu- 
tion for the 10° M GC simulations (C and F), which are shown in 
Fig. [4] For the 10 6 Mq GC simulations, we address the analysis of 
the two assumed specific mass-loss rates separately. Fig. [6] shows 
simulations A, D and G, which employ the higher mass-loss rate of 
a — 4 x 10~ 19 s -1 . Fig.plshows simulations B and E using the 
lower rate of a = 4 x 10~ s _1 . We finally compare the relative 
gas content for all simulations in Fig. [8] which shows the ICM mass 
as a fraction of the total injected gas with time. 



3.2.1 "Typical" mass GC simulations 

The 10 5 Mq GC simulations (C and F) represent a "typical" GC 
mass (Gnedin & Ostriker|1997| > and so the results are significant 
for the majority of the Galactic GCs we observe. 

Here, we investigate how the Halo density affects the ICM gas 
evolution. The ICM content within these clusters is strongly influ- 
enced by the density of the Galactic halo medium, with the final 
ICM mass in simulations C and F (Fig. |4| differing by nearly a 
factor of 6. Simulation C reaches a steady state by about 10 Myr 
and contains (i.e. material within r t ) about 1.3 Mq of gas for the 
remainder of the simulation. The ICM mass in simulation F starts 
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Figure 4. The ICM content over time for two 10 5 Mq GCs with spe- 
cific mass-loss rate a = 4 X 10~ 19 s — 1 , moving through a Galac- 
tic halo medium at 200 kms" 1 . Simulation C uses a Halo density of 
10 _26 gcm — 3 and F uses 10 — 27 gem -3 . The black dotted line is the 
injected gas mass against time. 



to reach a plateau at ~ 24 Myr then, after 4 Myr, continues to in- 
crease roughly linearly. The rate at which the gas mass increases 
is 4% of the total stellar mass-loss rate. At 10 Myr, F contains 
~ 5 M Q of ICM gas, which is 40% of the injected material, i.e. 
60% has been stripped away by this point. For C, 90% of the in- 
jected material has been removed by 10 Myr. By 50 Myr the total 
ICM content in F is 7.7 M (12.3% of injected material) whilst for 
C, an impressive 98% of the injected gas has been stripped from the 
cluster. 

Because efforts to detect gas within GCs focus primarily on 
hydrogen, we determine the H content of the ICM by multiply- 
ing the ICM mass by the hydrogen mass fraction (X = 0.7). At 
10 Myr, the H mass content in simulations C and F are 0.9 Mq 
and 3.5 Mq, respectively, whilst at 50 Myr, F contains 5.4 Mq of 
hydrogen (C still contains ~ 0.9 Mq). However, this is just the gas 
that resides within the cluster tidal radius. Fig.[5]shows that the ICM 
gas predominantly resides in a trailing tail of stripped material. 
Only a fraction of this total hydrogen gas resides within the core 
region. Selecting an arbitrary radius of 2 pc, the fraction of ICM 
gas within this region is 15% (0.13 Mq) for C and 9% (0.32 M Q ) 
forF. 

In addition to the morphological structure of the ICM, Fig. [5] 
also shows the ICM temperature to be lower than 10 4 K, indicating 
that hydrogen will be almost entirely neutral; indeed, using equa- 
tion (1 1) of Faulkner & Freeman ljT977j shows this to be the case. 
This result suggests that HI, rather than H a observations will pro- 
vide more insightful information about the ICM content in GCs. 

In summary, the predicted gas content of "typical" GCs is 
around the upper limit values given in the literature. However, due 
to the morphology of the ICM gas, the actual mass of gas that one 
would determine from observational data will clearly depend on the 
GCs direction of motion relative to us. 
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Figure 5. Simulation F: A cutaway of a 3D rendered image showing 
log (p/ g cm -3 ) (upper panel) and log (T/ K) (lower panel) of the ICM 
gas from a 10 5 Mg GC simulation at 10 Myr. These data are taken from 
the 3 rd grid that has a side length of 23.9 pc. 



In conclusion, ICM gas is continuously stripped from typical 
clusters and the residual gas content will prove particularly chal- 
lenging to detect with current radio telescopes. The bow shock tem- 
perature of around 10 6 K (or 0.08 keV), is consistent with the X- 
ray emission temperatures observed by Okada e t al.|p007) . 



3.2.2 High mass GC simulations 

The 10 6 Mq clusters represent the upper tail of the Galactic GC 
mass distribution. Removing gas from them poses a particular chal- 
lenge, since these clusters have much deeper gravitational poten- 
tial wells as well as more stars contributing to the ICM. Unlike 
the 10 5 Mq GC simulations, where the ICM content eventually re- 
mains constant with time, the 10 6 Mq simulations (A, B, D, E and 
G) show a continuous increase (Fig.|6]and Fig.|7](. This increase oc- 
curs almost linearly, at some fraction of the rate of injected stellar 
gas. It is clear that this linearity is brought about by a balance be- 
tween the force due to ram pressure and that due to gravity; where 
the rate of ram-pressure stripping is less than the injection rate of 




5 10 
time (Myr) 

Figure 6. The ICM content over time for three 10 6 Mq GCs with spe- 
cific mass-loss rate a = 4 X 10 -19 s — 1 , moving through a Galactic halo 
medium. Simulations A and D move at 200 km s -1 through Halo den- 
sities of 10 — 26 gem -3 and 10 _27 gcm -3 , respectively. Simulations D 
and G move though a Halo density of 10 — 27 gem - 3 with velocities of 
200 km s _1 and 300 km s — 1 , respectively. The black dotted line is the in- 
jected gas mass against time. 
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Figure 7. The ICM content over time for two 10 6 Mq GCs with spe- 
cific mass-loss rate a = 4 X 10 -20 s — 1 , moving through a Galactic halo 
medium at 200 kms -1 . Simulation B uses a density of 10 — 26 gcm~ 3 
and E uses 10 — 27 gem -3 . The black dotted line is the injected gas mass 
against time. 
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Figure 8. The ICM mass within different GC simulations as a fraction of 
the cumulative injected gas mass against time. 



the stellar material. This is demonstrated in Fig.[8] where the curves 
flatten out once this equilibrium is established. This happens on a 
timescale of about 5 Myr. The invariability of this equipoise over 
periods of Myr is nicely illustrated by simulations D, E and G, 
where the curves remain almost parallel with the injected mass of 
gas. For A and B, the ICM mass relative to the total injected stellar 
material begins to rise again after a while. This is simply due to 
the accumulation of the denser upstream Halo medium. All of the 
10 6 Mq simulations establish a bound ICM, where mass increases 
over time above the observed upper limits. At 10 Myr, simulation 
A contains 78.9 Mq, D, 107.1 M and G contains 97.9 M , or- 
ders of magnitude above observational limits. Simulations B and 
E, which adopt the lower value of a, contain 4.6 M© and 7.5 Mq, 
respectively, by 10 Myr. 

It is clear that even under extreme conditions (of Galactic halo 
density or GC velocity), the environment alone cannot prevent a 
build-up of bound ICM. The factor that produced the greatest in- 
fluence on the ICM content is the specific mass-loss rate. Adopting 
a lower value of a, affects the ICM evolution in two ways: firstly, 
it reduces the amount of gas being injected per unit time, trivially 
resulting in less gas to be removed. Secondly, for a given ram pres- 
sure, a lighter ICM will have a larger fraction lifted out of the GC 
potential. This latter effect is evident in Fig.[8]from the relative lo- 
cations of D and E, as well as A and B, where simulations D and 
E have the same ram pressure but different values of a. The same 
is true for A and B, which experience a higher ram pressure than 
D and E. The lower a values constitute a 30 — 40% reduction in 
the relative amount of remaining ICM per unit time. Considering 
the impact of one's choice of a and the restrictions on varying the 
environmental parameters, we are led to conclude that modelling 
the GC correctly and how mass is lost from the stellar population 
is a key factor in furthering our understanding of the evolution of 
ICM gas within GCs. 

Fig. [£] shows that the ICM gas content in 10 5 M GCs is af- 
fected by the Halo density to a greater degree than in the 10 6 Mq 



GCs (simulations A and D). In Fig. [8] we note that the initial mass 
of ICM gas compared to the injected stellar mass increases above 
1 (by up to 10%). This is due to the congestion of upstream Halo 
gas interacting with the injected stellar gas at early epochs. This 
is evident from the location of the peak of this initial increase, 
which occurs when we expect stripped central Halo-ICM gas to 
begin crossing r t (~ 0.1 Myr). This also implies that the ICM gas 
within r t for the 10 5 Mq GCs is refreshed with new material on 
this timescale. 



4 MULTI-MASS KING MODEL 

We move to a multi-mass model of the GC in order to better de- 
scribe both the distribution of mass loss as well as the prescrip- 
tion of the mass loss itself within our GC model. To date, the pre- 
scribed mass loss has been treated as essentially a free parame- 
ter, q, where values are guided from observed cluster properties 
( |Knapp et al.|1973[|TayIer & Wood|i975) . This blanket mass-loss 
rate is applied over the entirety of the GC. This approach is not 
entirely accurate, as it is the RGB and AGB stars that almost ex- 
clusively provide material for the ICM gas. Furthermore, mass loss 
in these stars, although still not well understood, has been stud- 
ied extensively. Hence, we can prescribe mass-loss rates according 
to relevant empirical formulae (Reimers 1975 , Vassiliadis & Wood 
|1993| respectively) that are derived from stellar properties, not the 
heterogeneous properties of a GC. The initial investigations pre- 
sented in Sections [2] and [3] use the King model, which is derived 
using the assumption that all the stars have the same mass. How- 
ever, a real GC consists of a stellar population spanning a range 
of masses, where radial distributions are subject to mass segrega- 
tion via two-body relaxation. This means that massive giant branch 
stars are concentrated towards the core compared to an average dis- 
tribution that is dominated by the more numerous lighter stars. In 
order to model the AGB and RGB stars as separate populations 
with different radial distributions and independent mass-loss rates, 
we require a multi-mass King model. A multi-mass King model is 
one in which the total mass and overall stellar density distribution 
is constructed out of several single-mass King models, each one 
made from a different stellar mass. 



4.1 Implementation of the multi-mass King model 

In order to construct a multi-mass GC model, one requires a stellar 
population created from an initial mass function (IMF). We use the 
Kroupa IMF l |Kroupa|200"2| >: Equation {XT} 

£0.019 23 form>O.5M , 

£0.038 (-^J form < 0.5 M© , 



dn 
dm 



where 

£ = 



(ID 



30.55 x 10 5 for the 10 5 M© GC , 
30.55 x 10 6 for the 10 6 M© GC , 



to obtain the number of stars and the total stellar mass in each pop- 
ulation bin. From this, we obtain a median stellar mass for that bin. 
We construct a multi-mass King model following the description in 
|Da Costa & Freeman| < |1976| > and |Capuzzo Dolcetta et al.| < |2005| >. 
In summary, this method first constructs a King model GC and fills 
the resultant gravitational potential with a number of so-called sub- 
King models, each representing a different mass bin in the stellar 



10 W. R. Priestley, M. Ruffert, M. Salaris 



o 



-5 



-10 



I 1 1 1 1 1 1 1 1 


1 i i i i 1 i i i i 1 i i i 


r 






- 
- 




- 
- 








King model 


population 4 \ 




_ — multi-mass 


- population 6 




_ — population 1 


- population 8 




_ — population 2 
i i i i 1 i i i i 


- population 10 
1 i i i i 1 i i i i 1 i i i 





-2-10 1 
log(r/r c ) 

Figure 9. The stellar density of a King GC and multi-mass GC model. The 
values are given as a fraction of the central density, p , of the King model. 
The dotted line is the density profile of the King GC profile and the solid 
black line is the sum of the component sub-King models, a sample of which 
are represented by the coloured lines. The AGB and RGB populations are 
denoted by populations 1 and 2, respectively. 
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Figure 10. A Hertzsprung-Russel diagram showing the zero age horizontal 
branch position of two stars that have experienced Reimers mass loss on the 
RGB with rj = 0.2 (red line) and r) = 0.4 (blue line). The dashed black 
lines outline the instability strip and the solid black line denotes the location 
of the RGB for both stars. 



population. These sub-King models are related to each other by as- 
suming an equipartition of kinetic energy between them: 



2 

miff] 



mitJi where i > 2 



(12) 



We arbitrarily split our stellar population into ten mass bins where 
the heaviest two represent the AGB and RGB stellar populations 
(bins 1 and 2, respectively). The rest are split among the sub- 
giant branch and main sequence stars. A comparison of the King 
and multi-mass GC models is shown in Fig. [9] Using a theoretical 
isochrone from the BaSTI database ( Pietrinferni et al. 2004), we 
normalise the IMF such that the resulting stellar population repre- 
sents a currently evolving GC. In order to remain consistent with 
the gas physics used in VF77, we employ an isochrone with scaled 
solar metallicity (Z — 0.002) at an assumed age of 13 Gyr. This 
yields a GC where stars more massive than the AGB population 
(~ 0.83 M ) no longer exist within the cluster (except for stellar 
remnants such as neutron stars and white dwarfs). We also use this 
isochrone to determine the stellar properties of this population. 

As stated before, the purpose of using a multi-mass King 
model is to better represent the mass loss that occurs within GCs. 
Each of the RGB and AGB stellar populations therefore have their 
own specific mass-loss rate, qrgb and ciagb. To achieve this, 
we appropriately sample the isochrone and obtain random stel- 
lar properties for each star in our RGB and AGB population bins 
(Equation These are then inserted into the relevant empirical 
mass-loss formula (Reimers 1975, Vassiliadis & Wood 19931 to 
obtain individual mass-loss rates. It is important to note that the 
Reimers mass-loss formula has a free parameter, -q, which typi- 
cally adopts the values 0.2 or 0.4, where the larger value leads 
to a higher rate of mass loss. The RGB and AGB specific mass- 
loss rates are determined by dividing each population's combined 



mass-loss rate by its total stellar mass. The equivalent global a 
value is determined by dividing the overall mass-loss rate (i.e. from 
both the AGB and RGB populations) by the GC mass. It is in- 
teresting to ask how such an a value compares to that used by 
VF77. With r] = 0.2, « = 4 x 10~ 20 s~\ an order of magni- 
tude less than that adopted by VF77 and coincidentally matching 
the observationally determined value in jTayler & Wood| ( |1975| >. 
rj = 0.4 yields a global a of roughly 1 x 10 _iy s _i , one fourth 
the value adopted by VF77. Adopting rj = 0.2 leads to 0.1 Mq of 
gas lost during the RGB phase and a total cluster mass-loss rate of 
1.3 x 10~ 6 M yr _1 (for a 10 6 M GC), of which 58% is from 
the RGB population and 42% from the AGB. Taking rj = 0.4 
leads to 0.2 M lost on the RGB and a total mass-loss rate of 
3.5 x 10" 6 M yr -1 , 85% of which is due to t he RGB stars and 
15% is from the AGB. |McDonald et alT| l |2009^ determine a total 
cluster mass- loss rate of 1.2^;° x 1O~ 5 M yr" 1 for uj Centauri 
and suggest that 0.2 — 0.25 M of gas is lost on the R GB. Using 
a cluster mass of 3.1 ± 0.3 x 10 6 M |Miocchi|20"Tb] , this corre- 
sponds to con = 1.2lo'g x 10 _19 s~ , which favours higher 
mass-loss rates, but is consistent with our derived a values. 

The resultant zero age horizontal branch (ZAHB) position of 
a star having undergone our two different mass-loss rates is shown 
in Fig. [10] The diagram shows how rj = 0.4 produces a blue ZAHB 
star with very little of its envelope remaining. This implies that the 
value of a adopted by VF77 is perhaps too extreme, resulting in a 
ZAHB star that will have lost all of its envelope. An a derived from 
observations is cluster dependent and requires an adopted mass-to- 
light ratio; this results in an unrelated range of possible values that 
are hard to apply generally (c.f. the values in Knap p et al.|1973) 
|Tayler iSTWo od 1975 ). We prefer the approach adopted here; it is 
better suited for modelling mass loss in GCs and investigating the 
subsequent ICM evolution in a consistent and methodical way. 
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Table 5. Summary of the simulations that compare the ICM evolution using 
a King model GC and a multi-mass GC. 
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Figure 11. The radial profile of the overall local mass-loss rate for the multi- 
mass (solid line) and King (dotted line) GC models. The mass-loss rates are 
normalised to the King model's central mass-loss rate. 



With consistently determined specific mass-loss rates, we can 
test how modifying the distribution of mass loss (i.e. the difference 
between using a King model and a multi-mass GC model) affects 
the GC ICM evolution. For the multi-mass model, the RGB and 
AGB stars are more centrally concentrated, resulting in a larger 
proportion of gas being injected deeper within the potential well of 
the cluster, as shown in Fig.[TT] Therefore, for a given ram pressure, 
we expect less gas to be removed from the GC potential compared 
to the King-model simulation. The result, is a poorer agreement 
between our simulations and observations! As the first-ever inves- 
tigation describing mass loss both in terms of empirical laws and 
distribution within the cluster, it is essential to ascertain how much 
the ICM evolution depends upon the GC model. The subsequent 
evolution of these two cluster models moving through a Galactic 
halo medium is discussed in Sectionl4~2l 



4.2 RGB and AGB mass-loss simulations 

We run two high-resolution simulations each for a 10 5 Mq (64 3 
with seven nested grids) and 10 6 Mq (64 3 with six nested grids) 
GC in order to compare the difference between using a King model 
GC and a multi-mass model GC (i.e. four simulations in total); 
these are summarised in Table [5] The method used to produce the 
initial King model that is involved in constructing our multi-mass 
models, is slightly different to the method employed for the King 
models described in Sections [2] and [3] Therefore, the King models 
used in this section and previous sections are comparable, but not 
identical. The parameters that describe the new King models are 
provided in Table [6] We run a new set of King-model GC simula- 
tions in order to maintain a consistent comparison with our multi- 
mass GC simulations. The initial conditions for all simulations are 
the same as in models D and F (10 6 Mq and 10 5 M GCs respec- 
tively). That is, a Halo temperature of 10 K, a Halo density of 



Simulation 
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GC model 


a (s" 1 ) 
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multi-mass 
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10 s 


King 
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Table 6. Properties of the two King models created for the purpose of con- 
structing the multi-mass GC models. 
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1.03 x 10 5 


1.653 x 10 


-18 


0.51 


24.5 


1.03 x 10 6 


1.653 x 10 


-17 


0.51 


24.5 



tively easy to strip gas from the 10 5 Mq GCs, we use the Reimers 
mass-loss law with r\ — 0.4 (corresponding to a — 1 x 10~ 19 s _1 ). 
For the 10 6 Mq GC, we use the lower extreme of t] = 0.2 
(a — 4 x 10 _20 s _1 )in order to maximise gas stripping. In the 
following discussion, we compare how the ICM evolution changes 
between the King GC model and the multi-mass description of the 
cluster. We first compare the 10 5 M GC simulations (J and K, 
Fig. 1 1 2^ and then discuss the results from the 10 6 Mq simulations 
(H and I, Fig,|13[>, Lastly, we compare the relative evolution of the 
ICM mass for all sim ulations (Fig.|14|l. 



In Section 



3.2 



the long-term build-up of ICM gas in 10° M Q 
clusters was prevented, due to the ram pressure experienced as the 
GC moved through the ambient medium. Fig.[l2]shows this is also 
the case for the multi-mass cluster. Compared to the King model 
(J), however, more gas is retained by the GC potential before the 
ICM content plateau (at ~ 24 Myr). The plateau mass of gas within 
the cluster for the rest of the simulation is roughly 1.96 Mq, fluc- 
tuating by up to 9%. This fluctuation is due to the central reservoir 
of gas increasing in size until Kelvin-Helmholtz (KH) instabilities 
develop on the surface tangential to the flow of Halo gas around the 
ICM gas. These KH instabilities result in an episode of enhanced 
gas stripping, which reduces the size of the core. This cycle of ICM 
growth and shedding occurs on timescales of 13 ± 1.5 Myr. The 
King model reaches a plateau of 0.97 Mq in half the time and can 
vary by up to 6%. In this steady-state condition, the multi-mass 
simulation has just over twice as much gas compared to the King- 
model simulation. The corresponding hydrogen content within r t 
for these two model clusters are 0.68 Mq and 1 .37 M for J and K, 
respectively, which are above observational limits. However, within 
a core region of 2 pc, there is as little as 0.15 Mq and 0.29 Mq of 
hydrogen gas present, about the same order as the most stringent 
observational upper limits, but much less than most upper limits. 
Again, most of the gas within r t , resides within the tail of stripped 
material, where photometric detection will be line-of-sight depen- 
dent. 

For the 10 6 Mq model simulations, the multi-mass GC also 
retains more gas than the King-model simulation. Fig. [T3 
that simulation H develops similarly to simulation E (Fig 



shows 
7), the 

rate of gas stripping due to ram pressure being stable over tens 
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Figure 15. A cutaway of 3D rendered images showing log (p/ gem -3 ) (upper panels) and log (T/ K) (lower panels) of the ICM gas at lOMyr for 
10 5 Mq GC represented by a King model (simulation J; left panels) and multi-mass model (simulation K; right panels). These data are taken from the 3 r 
grid that has a side length of 24.5 pc. 




Figure 16. A cutaway of 3D rendered images showing log (p/ g cm -3 ) (upper panels) and log (T/ K) (lower panels) of the ICM gas at 10 Myr for 
10 6 GC represented by a King model (simulation H; left panels) and multi-mass model (simulation I; right panels). These data are taken from the 3 r 
grid that has a side length of 24.5 pc. 
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Figure 12. The ICM content over time for two 10 5 M GCs moving 
through a Galactic halo medium, one modelled with a King model (J) and 
the other with a multi-mass model (K). The representative stellar popula- 
tions have the same total mass-loss rate (equivalent to a specific mass-loss 
value of a = 1 X 10~ 19 s _1 ). The black dotted line is the injected gas 
mass against time. 
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Figure 13. The ICM content over time for two 10 6 M GCs moving 
through a Galactic halo medium, one modelled with a King model (H) and 
the other with a multi-mass model (I). The representative stellar populations 
have the same total mass-loss rate (equivalent to a specific mass-loss value 
of a = 4 X 10 — 20 s -1 ). The black dotted line is the injected gas mass 
against time. 



of Myr. This results in the ICM mass increasing linearly at about 
57% the rate of injection of stellar material. The multi-mass model 
(simulation I) experiences a much lower rate of gas stripping, with 
the ICM increasing at about 84% of the rate of stellar mass loss 
(47% more than H). Since the ICM content increases continu- 
ally with time, we again choose to assess the hydrogen content at 
10 Myr. Model H contains 5.0 M Q within r t , of which 4.6 M re- 
sides within the central 2 pc. In comparison, simulation I contains 
7.7 Mq of hydrogen, 7.6 M residing within 2pc. Therefore, we 
expect to observe a central reservoir of gas for massive clusters 
moving through the tenuous upper Halo. 

The relative ICM content, compared to the injected mass, from 
all simulations is shown in Fig. [14] It highlights well, the signifi- 
cant influence that cluster mass has on the mass of retained ICM. 
From these simulations, it appears that removing gas from the grav- 
itational potential of massive GCs may be a bigger challenge than 
previous studies have inferred. The ~ 50% difference (100% in the 
case of 10 J Mq clusters) in ICM content generated by altering the 
distribution of stellar mass loss re-enforces the conclusion of Sec- 
tion[3] it is important to appropriately model the stellar population 
and associated mass loss. 

We present a snapshot at 10 Myr of the gas density and t ern- 
16 M 



15 



and Fig. 



16 



perature for the 10° M and 10° M GCs in Fig 
respectively. These figures show the morphological differences that 
emerge from using the King model (left panel of Fig. [15] and Fig. 
|16| l and the multi-mass one (right panel). The multi-mass GC sim- 
ulations, show a larger volume of denser material in the core com- 
pared to the King-model simulation. The 10 5 M simulations show 
a cool tail of stripped gas (which contains neutral hydrogen) with 
no evidence of a distinct, bound core of gas in the central regions. In 
contrast, the 10 6 M simulations show a central mass of gas that is 
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Figure 14. The ICM mass within different GC simulations as a fraction of 
the cumulative injected gas mass against time. 
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thermally distinct from the tail of stripped material, indicating that 
this gas is bound and cools over time. The hydrogen in this dense 
central reservoir is atomic. 



5 DISCRETE STELLAR POPULATION CLUSTER 
MODEL 

Section |4~2] highlighted the influence the GC model has on simu- 
lating ICM evolution. The King and multi-mass GC models both 
represent the stellar population as a continuous stellar density dis- 
tribution. This approximation is upheld in the limit of a large num- 
ber of stars per unit volume. This is indeed the case (although is 
less valid at the cluster outskirts) for the King model, which has 
a homogeneous stellar population. For a multi-mass model, this 
approximation remains valid for the numerous (10 5 — 10 6 ) main- 
sequence stars. However, the approximation of a continuous stel- 
lar distribution breaks down for the RGB and AGB populations, 
since they account for less than a percent of the total population. 
For our 10 5 M Q GC, the IMF produces 814 RGB stars and just 
2 AGB stars. The 10 6 M GC contains 8140 and 20 of each, re- 
spectively. Modelling these populations as individual stars orbiting 
within the GC potential will alter the distribution of gas injection 
from a monolithic, symmetric one to a time-dependent distribution 
of asymmetrically dispersed points. In such a distribution, the ICM 
injection is not centred on the deepest parts of the gravitational po- 
tential, allowing a larger fraction to be liberated by ram pressure. 
A distribution of independent sources of material, increases the ef- 
fective surface area upon which ram pressure will act to lift gas 
out of the potential. Furthermore, an asymmetric ICM morphology 
will provide more seed perturbations from which KH instabilities 
can grow, leading to the mixing and fragmentation of the ICM. The 
importance of KH instabilities on the fragmentation and break-up 
of a cloud moving through an ambient medium, has already been 
demonstrated in Agertz et al. (2007 1. Although our ICM gas is con- 
tinuously replenished, we hypothesise that these effects will have a 
significant impact on ICM evolution, leading to a greater level of 
gas stripping. In the following sections, we present our discrete GC 
model and subsequent simulations. 

5.1 Rationale, setup and initial conditions 

For the purpose of numerical efficiency, we make a few simplifying 
approximations in modelling our discrete GCs. We choose to repre- 
sent our RGB and AGB populations as point sources of mass loss, 
moving as test particles within the smooth GC potential. Since we 
are only concerned with the stellar winds of these two populations, 
modelling the other stars as discrete objects is unnecessary. Sec- 
ondly, by neglecting full N-body calculations, we avoid consuming 
CPU time on resolving unrelated binary interactions. Furthermore, 
ignoring star-star interactions, isolates the influence of the discreti- 
sation and stellar motion (i.e. we don't include the effect of a noisy, 
time- varying potential). 

We build our model by giving each star, from each population, 
an initial position and velocity within our simulation volume. We 
use the multi-mass GC model in spherical polar coordinates as the 
framework for obtaining these initial values via Monte-Carlo tech- 
niques. The angular positions of our stars are distributed evenly, 
whereas the radial position is weighted according to the sub-King 
model for each population. At each radial position, the magnitude 
of the stellar velocity is taken to be the local root mean squared 
velocity and the vector direction is chosen at random. Each star is 



given its own individual empirical mass-loss rate (the same as that 
used in calculating the specific mass-loss rates in Section]?}. 

The procedure used to employ the discrete-model GCs into 
the simulations is slightly modified compared to that described in 
Section[2] For the injection of stellar material, each cell no longer 
holds a stellar density value. Instead, we track the position of each 
star, relative to the grid and add gas directly into the cell it resides 
in. Cells that contain no stars do not have material injected into 
them. The local ICM mass, momentum and energy properties are 
updated accordingly (Section|2j. The local momentum and energy 
conservation Equations, Q and j5} respectively, must be gener- 
alised slightly to include the extra terms that vanish in the limit of 
an isotropic, continuous stellar distribution. The local momentum 
is conserved according to 

Pn + lVn+l = PinjV* + p n V n , (13) 

and the energy injected by the stellar wind, U a , is modified such 
that 

1 2 1 2 

U a = p+ -v* - V*.V n+ l + 2 U n+i • (14) 

After injecting mass onto the grid, the stars are moved (with the 
hydrodynamic timestep) under the influence of gravity, using the 
Euler method and a predictor-corrector scheme. The rest of the hy- 
drodynamic calculations follow the sequence described in Section 

m 

5.2 Evolution of the ICM within a discrete stellar 
environment 

For each GC mass, we run a total of five low-resolution simula- 
tions using three nested grids, L to P for the 10 6 Mq GC and Q 
to U for the 10 5 Mq GC. For all simulations, the initial conditions 
for the gas on the grid are p H = 10~ 27 gcm -3 , T H = 10 5 ' 5 K 
and vac = 200 km s -1 . Each simulation has a different set of 
randomly determined initial positions and velocities for the stars. 
The Halo initial conditions remain the same as those used in previ- 
ous sections. We run five simulations in order to see what effect, if 
any, the choice of initial stellar positions and velocities has on the 
global ICM evolution. In this section, we present the results from 
the 10 5 Mq simulations, followed by a brief description of those 
from the 10 6 Mq simulations. We then discuss some of the techni- 
cal challenges faced when attempting high-resolution simulations 
and conclude with a summary of the main results from this experi- 
ment and discuss the consequences for future investigations. 

Section |4~2| showed that a multi-mass GC model (employing 
mass-loss laws) results in a higher retention of ICM gas than the 
King model. The predicted ICM mass from our discrete multi-mass 
GC simulations are given in Fig. [T7] Four of the five discrete GC 
simulations show a plateau mass of 0.92 Mq, varying by about 
10% due to stellar orbits moving gas injection around the potential 
well. This is 60% less than a low-resolution continuous multi-mass 
GC model (2.28 M Q ± 12.5%) and 35% less than an equivalent 
King-model simulation (1.43 Mq ± 1%). The discrete model that 
predicts the highest ICM mass (T) is comparable with the King- 
model simulations at 1.35 Mq ± 10%. This simulation retains more 
gas from the start, implying that some members of its stellar pop- 
ulation are in lower eccentricity orbits, spending more time deeper 
in the potential well compared to other simulations. However, stel- 
lar orbits have not been traced throughout the simulations and this 
hypothesis has not been confirmed. 

The log density and temperature profiles of these simulations 
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Figure 18. A cutaway of 3D rendered images showing log (p/ gem -3 ) of the ICM gas at lOMyr for a 10 5 Mg GC, represented by discrete multi-mass 
GC models that have different initial stellar positions and velocities (simulations Q, R, S, T, and U). The bottom-right panel is from an equivalent continuous 
multi-mass simulation. These data are taken from the 3 rd grid which has a length of 24.5 pc. 



at 10 Myr are presented in Fig.|18|and Fig.[T9] respectively. Com- 
pared to a low-resolution continuous multi-mass GC simulation 
(bottom-right panel in Fig. [T8] and Fig. 1 19[>, the ICM morphology 
from the discrete simulations show less ICM content due to a nar- 
rower cross-sectional profile. The streaming tail of gas shows more 
structure and a lower density. A few stars are also observed losing 
mass separately from the main region of ICM gas. A summary of 
the amount of hydrogen within the central 0.2 pc is given in Ta- 
ble [7] These values range from 0.11 Mq to 0.17 Mq and reside 
between the 0.1 Mq value for the King-model simulations and the 
0.29 Mq for the continuous multi-mass model. 

A discrete model for the 10 6 Mq GC simulations produced no 
change in the overall ICM mass content with time, compared to a 
low-resolution continuous multi-mass simulation. This is likely be- 
cause the stellar population is unresolved in these low-resolution 
simulations. Accordingly, the discrete simulations behave much 
like the continuous stellar density case. However, it is also possi- 
ble that, for the massive GC, the number density of RGB and AGB 



Table 7. The mass of hydrogen residing within the core of the 10 5 Mq GC 
discrete multi-mass simulations. 



Simulation H mass within 2 pc ( Mq) 



Q 0.11 

R 0.12 

S 0.15 

T 0.17 

U 0.16 



stars is high enough to warrant the use of a continuous stellar den- 
sity distribution. To determine the effect of discretising the stellar 
population on the ICM evolution, the stellar population needs to be 
resolved. This requires higher resolution simulations with a larger 
number of nested grids so that the individual stars can be sepa- 



simulation R 
discrete multi-mass 




ios(r/K) 
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Figure 19. A cutaway of 3D rendered images showing log (T/ K) of the ICM gas at 10 Myr for a 10 5 Mq GC, represented by discrete multi-mass GC 
models that have different initial stellar positions and velocities (simulations Q, R, S, T, and U). The bottom-right panel is from an equivalent continuous 
multi-mass simulation. These data are taken from the 3 rd grid which has a length of 24.5 pc. 



rated on the grid. Nevertheless, it is encouraging that these initial 
10 5 Mq GC simulations indicate that discretising the stellar pop- 
ulation has a favourable effect on reducing the level of ICM gas 
within GCs. 

We attempted high-resolution simulations (64 3 with six lev- 
els of refinement). However, we discovered a numerical problem in 
our simple discrete GC model. In these simulations extreme tem- 
peratures spanning ~ 10 9 K to ~ 0.1 K are observed. We realised 
that by injecting the gas into a single cell, we create a numeri- 
cally unresolved injection process and do not correctly simulate 
the changes in the gas flow properties as the stellar wind inter- 
acts with the Halo gas. This discontinuous peak in density, con- 
tained within a solitary cell, is an obstacle that Halo gas cannot 
flow around properly. This leads to artificially low densities in the 
star's wake (p ~ 10 -30 gem -3 ) as gas leaves these cells due to 
the relative motion of the Halo. This wake gas is not replenished by 
the upstream wind material at an appropriate rate since the flow of 
this material is not established over enough cells. This wake gas not 



only has a reduced density but is also artificially heated to ~ 10 K, 
causing numerical errors later on in the simulation such as ex- 
tremely cold gas. This cold gas has a low internal energy which 
results in tiny energy injection timesteps and the simulation grinds 
to a halt. This numerical effect, caused by injecting gas into a single 
cell, becomes more extreme with increasing resolution and so our 
initial simulations therefore do not show temperature variations be- 
yond that expected from cooling. Although the gas physics in these 
low-resolution discrete simulations suffer from numerical artefacts, 
we firmly believe that this does not change the predicted low lev- 
els of ICM gas. The continuous stellar density simulations do not 
suffer from this numerical issue because the gas is injected over the 
entire cluster and the flow of this new gas is extended over most of 
the cells in the simulation volume. Future work will present sim- 
ulations where the injection of stellar wind material from discrete 
sources is spread over several cells. The resultant high-resolution 
simulations aim to confirm our initial findings, as well as determine 
the effect that discretisation has on the ICM evolution of 10 6 Mq 
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Figure 17. The ICM content over time for seven 10 5 Mq GCs moving 
through a Galactic halo medium. Five are represented by a discrete multi- 
mass model with different initial stellar positions and velocities (simula- 
tions Q through U; solid coloured lines), one is represented by a King model 
and another is represented by a continuous multi-mass model. The represen- 
tative stellar populations in each GC model have the same total mass-loss 
rate (equivalent to a specific mass- loss value of a = 1 X 10 -19 s" 1 ). The 
black dotted line is the injected gas mass against time. 

GCs. 

Our initial findings from using a discrete GC model, are pre- 
sented here as a proof-of-concept result. Breaking up the injection 
of stellar winds into localised sources, enables gas to escape the 
cluster potential more easily, as well as encourage the growth and 
development of instabilities that lead to the fragmentation and dis- 
sipation of a replenished intra-cluster medium. 



6 SUMMARY AND DISCUSSION 
6.1 Discussion 

Our work highlights the importance of carefully modelling the 
amount and distribution of mass loss, when simulating the evo- 
lution of ICM gas in GCs. Defining the nature of mass loss from 
population II giants is therefore essential for understanding the evo- 
lution of gas within GCs. In this work, we employed the Reimers 
formula, which is often applied in many fields to describe mass 
loss from GC RGB stars. However, the Reimers law is calibrated 
on population I stars and isn't successful at explaining all the ob- 
servations of GC stellar populations. There are currently only a 
handful of alternative formulae available that have been calibrated 
against globular cluster RGB and AGB stars. They give quite differ- 
ent mass-loss rates, some suggesting that these winds are episodic 
in nature, highlighting the fact that the mechanisms governing 
mass loss are poorly understood (Origlia et al. 2007; Meszaros 
|etal.|20 09; McDonal d et al.|2009| >. Although on average the same 
amount of gas must be injected into the ICM, episodic mass-loss 
may have an encouraging role in removing gas from the cluster po- 



tential. Episodic mass-loss will cause variations in the local stellar 
mass-loss rate due to some stars turning off their mass loss whilst 
others begin their duty cycle. In those regions where more stars are 
in an "off" phase, a larger portion of ICM gas could be removed 
by ram pressure per unit time. However, simulation is the only way 
to determine the influence this sort of mass loss will have on the 
evolution of the GC ICM. 

For "typical" (i.e. 10 Mq) GCs, we used the higher value 
for the Reimers stellar mass-loss rate (i.e. rj = 0.4) as described 
in Section [4~2] Even with this higher rate, our predicted ICM con- 
tent is less than the majority of observational upper limits. There- 
fore, if we adopt a lower Reimers mass-loss rate {q = 0.2) for 
stars on the RGB, we would obtain an even lower predicted ICM 
mass. In contrast, the massive (10 6 Mq) GC models used the lower 
Reimers mass-loss rate and yet still show a disparity between the 
predicted and observed ICM mass for the massive GC models. We 
propose that discretising the distribution of mass loss in simula- 
tions will improve the gap between theory and observation. How- 
ever, for 10 6 M© clusters employing higher stellar mass-loss rates, 
it is likely that other mechanisms, such as those discussed in Sec- 
tion [T| will need to be included in order to prevent the build-up of 
ICM gas within massive GCs. 

We note that the radial distribution of our stellar population 
is fixed and does not change with time. Real GCs, however, are 
dynamical systems in quasi-equilibrium states, where stellar distri- 
butions change with time. Events such as the hardening of binary 
systems and stellar mass loss, result in an expansion of the core; 
this would make ram-pressure stripping more effective. Since we 
model mass loss, we must determine whether our adopted static 
(King and multi-mass) models are justified. To do this, we look 
at the stellar evolutionary and GC central relaxation timescales, as 
well as the fraction of stellar mass lost via winds between Disk 
crossing events. The timescales over which our RGB and AGB stars 
evolve are of the order 10 s years. Similarly, the central rel axation 
timescales of Galactic GCs show a median of 3.4 x 10 s yr |Harris| 
|1996[ >. The central relaxation timescale of our adopted King mod- 
els is roughly 9.0 x 10 7 yr (with a half-mass relaxation timescale 
of 2.6 x 10 8 yr). These numbers indicate that there may be small 
changes in the stellar distribution over the course of the simulations, 
implying that static models provide an upper limit on the predicted 
ICM content. In 10 8 yr, however, the total change in mass of the 
RGB and AGB populations due to stellar winds is only 2%, which 
is negligible. In our simulations, the effect of dynamical heating of 
the stellar distribution due to mass loss can therefore be ignored 
and the use of static models is appropriate. 

6.2 Summary 

We validate a 3D hydrodynamics simulation and undertake numer- 
ical investigations on the effect of the ICM interaction with the 
Galactic halo medium. We find that the GCs motion through the 
Halo, is sufficient in limiting the build-up of gas within "typical" 
mass globular clusters (i.e. 10 Mq), to levels that are below most 
observational limits. For massive GCs (10 6 Mq), however, we pre- 
dict the steady build-up of a cool reservoir of ICM that should be 
readily observable. We present simulations that employ a multi- 
mass GC model and incorporates empirical mass-loss formulae to 
model the RGB and AGB stellar winds. The use of these mass- 
loss formulae provides a way to consistently model the rate of gas 
injection into our GC ICM. The rates obtained from these formu- 
lae, show that previously assumed mass-loss rates were very high 
and result in blue HB stars. The implementation of a multi-mass 
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GC model, where the giant branch stars reside deeper in the grav- 
itational potential well, results in the retention of more ICM gas 
than previously predicted with King models. Finally, we present 
some preliminary work in which we model the GC as a discrete, 
orbiting stellar population. Although presented in its exploratory 
stages, our low-resolution 10 5 Mq GC simulations show that this 
approach readily predicts an ICM mass that is within observational 
limits. This result indicates that assuming a continuous stellar den- 
sity distribution may overestimate the predicted GC gas content by 
a significant fraction. We summarise our main results below: 

• moving GC through the Galactic halo medium reduces ICM 
content: 

- 10 5 Mq : ICM stripped to within observational limits 

- 10 6 Mq : ICM mass increases steadily with time 

• the mass of retained ICM gas is sensitive to the specific 
mass-loss (a) parameter used: 

- empirical mass-loss laws give lower a than that used in 
previous investigations; 

• more ICM is retained in a multi-mass GC model than in a King 
model GC; 

• discretising the stellar population amplifies and increases 
the influence of ram-pressure stripping. 
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